############## Step-by-step SCCs 
# rho/eta to use
rho = output.3pct$rho
eta = output.3pct$eta

##### Read in raw marginal damages of final results ----
# This is TxB (30x10k). 30 time steps: 2005:10:2295. 10k draws.
first.row.dice = 2 # keep 2015 on to match handling of MSW/UW merged data (which starts in 2021)

# Marginal damages
md.dice = as.matrix(fread(data.dir%&%'/DICE/part 5-6/DICE_part5-6_MC_MSW_IPF_PW05_2015_marginal_damages.csv'))
md.dice = md.dice[first.row.dice:nrow(md.dice),]

# Growth rates
md.g.dice = as.matrix(fread(data.dir%&%'/DICE/part 5-6/DICE_part5-6_MC_MSW_IPF_PW05_cpc_growth.csv'))
md.g.dice = md.g.dice[first.row.dice:nrow(md.g.dice),]
# convert from discrete rates to continuous rates

T.dice = nrow(md.dice)
B.dice = ncol(md.dice)

#  DICE is run on ten year timesteps from 2005:10:2295
timestep.dice = seq(2005, 2295, by=10)
timestep.dice = timestep.dice[first.row.dice:length(timestep.dice)]
n.years.dice = rep(10, length(timestep.dice))
length(timestep.dice)

rownames(md.dice) = timestep.dice
rownames(md.g.dice) = timestep.dice

### Replicate SCCs
st.yr.dice = 2015

# Calculate SCCs
calc.scc = function(rho=rho, eta=eta, md, g, timestep=10, all=FALSE, mean=TRUE, break.cor=FALSE) {
  r = rho+eta*g
  r[1,] = 0 # discount to period 1, so no discounting in that period
  
  DF = apply(r, 2, function(x) exp(-cumsum(x)*timestep))
  
  PV.MD = DF*md*timestep
  SCC = apply(PV.MD, 2, sum)
  E.SCC = mean(SCC)
  Med.SCC = median(SCC)
  
  if (break.cor) {
    if (mean) scc = sum(apply(timestep*md, 1, mean)*apply(DF, 1, mean)) # mean DF for each year times mean dmg for each year
    if (!mean) scc = sum(apply(timestep*md, 1, median)*apply(DF, 1, median)) # mean DF for each year times mean dmg for each year
    return(scc)
  }
  if (all) return(SCC)
  if (!all) {
    if (mean) scc = E.SCC
    if (!mean) scc = Med.SCC
    return(scc)
  }
}

# Inflation adjustment to 2020
# GDP deflator for inflation adjustments
infl = fread(data.dir%&%'/GDPDEF.csv')
setnames(infl, tolower(names(infl)))

mc.infl.adj = infl[year==2020,gdpdef]/infl[year==2011,gdpdef]

# Calculate SCC with Ramsey
sccs.ramsey = calc.scc(rho=rho, eta=eta, md=md.dice, g=md.g.dice, all=TRUE)*mc.infl.adj
mean(sccs.ramsey) # $78/ton

# Calculate SCC with Constant discounting
sccs.const = calc.scc(rho=0.03, eta=0, md=md.dice, g=md.g.dice, all=TRUE)*mc.infl.adj
mean(sccs.const) # $136/ton

# Calculate SCC with uncertain r and g, but breaking their correlation:
sccs.ramsey.break.cor = calc.scc(rho=rho, eta=eta, md=md.dice, g=md.g.dice, all=TRUE, break.cor=TRUE)*mc.infl.adj
sccs.ramsey.break.cor.med = calc.scc(rho=rho, eta=eta, md=md.dice, g=md.g.dice, all=TRUE, break.cor=TRUE, mean=FALSE)*mc.infl.adj
sccs.ramsey.break.cor # ~$10k without reflecting the correlation
sccs.ramsey.break.cor.med # $66

# Row 4: SCC with average MSW growth rate
md.dice.pt4 = as.matrix(fread(data.dir%&%'/DICE/part 4/DICE_part4_MC_MSW_IPF_PW05_2015_marginal_damages.csv'))
md.dice.pt4 = md.dice.pt4[first.row.dice:nrow(md.dice.pt4),]
rownames(md.dice.pt4) = timestep.dice

md.g.dice.pt4 = as.matrix(fread(data.dir%&%'/DICE/part 4/DICE_part4_MC_MSW_IPF_PW05_cpc_growth.csv'))
md.g.dice.pt4 = md.g.dice.pt4[first.row.dice:nrow(md.g.dice.pt4),]
rownames(md.g.dice.pt4) = timestep.dice

# Deterministic SCC with constant 3% (part 4)
sccs.const.pt4 = calc.scc(rho=0.03, eta=0, md=md.dice.pt4, g=md.g.dice.pt4, all=TRUE)*mc.infl.adj
median(sccs.const.pt4) # $56/ton
 

# Side calculation: Deterministic growth with Ramsey
sccs.ramsey.pt4 = calc.scc(rho=rho, eta=eta, md=md.dice.pt4, g=md.g.dice.pt4, all=TRUE)*mc.infl.adj
mean(sccs.ramsey.pt4) # $68/ton

#### Report SCCs for 6 scenarios, both mean and median ----
# Dimensions of SCC variation
# Scenarios (6)
# Mean/median (2)
scc = array(NA, dim=c(6,2), dimnames=list(c(1:6),c('Mean','Median')))

## Scenario 1: Actual IWG
# IWG was in 2007$, so we need a different inflation adjustment
iwg.infl.adj = infl[year==2020,gdpdef]/infl[year==2007,gdpdef]

## This is IWG result is for a 2015 pulse, yielding $34/ton
library(openxlsx)
scc.iwg.dice = read.xlsx(data.dir%&%'/DICE/part 1/EPA_DICE_2015_SCC_0.03.xlsx')

scc[1,'Mean'] = mean(as.numeric(as.matrix(scc.iwg.dice)))*iwg.infl.adj
scc[1,'Median'] = median(as.numeric(as.matrix(scc.iwg.dice)))*iwg.infl.adj

quantile(as.numeric(as.matrix(scc.iwg.dice))*iwg.infl.adj, probs=c(0.025, 0.975))

## Scenario 2: IWG + deterministic emissions
scc.dice.pt2 = as.matrix(fread(data.dir%&%'/DICE/part 2/DICE_part2_MC_2015_scc_values.csv'))

scc[2,'Mean'] = mean(scc.dice.pt2)*iwg.infl.adj
scc[2,'Median'] = median(scc.dice.pt2)*iwg.infl.adj

quantile(scc.dice.pt2*iwg.infl.adj, probs=c(0.025, 0.975))

## Scenario 3: IWG + deterministic emissions + average UW population growth
scc.dice.pt3 = as.matrix(fread(data.dir%&%'/DICE/part 3/DICE_part3_MC_2015_scc_values.csv'))

scc[3,'Mean'] = mean(scc.dice.pt3)*iwg.infl.adj
scc[3,'Median'] = median(scc.dice.pt3)*iwg.infl.adj

quantile(scc.dice.pt3*iwg.infl.adj, probs=c(0.025, 0.975))

## Scenario 4: IWG + deterministic emissions + average UW pop + average MWS growth
scc[4,'Mean'] = mean(sccs.const.pt4)
scc[4,'Median'] = median(sccs.const.pt4)

quantile(sccs.const.pt4, probs=c(0.025, 0.975))

# Scenario 4.5: Same but with rho/eta
# The result is quite similar if we use our discounting parameters 
# ($\rho=0.8\%$, $\eta=1.53$) before adding stochastic economic growth 
# from MSW--the mean and median SCC values are \$68 and \$59 respectively. 
# This illustrates that, without substantial uncertainty in economic growth, 
# the application of the discounting rule does little on its own.
mean(sccs.ramsey.pt4) # $68
median(sccs.ramsey.pt4) # $59

# Scenario 5: Stochastic growth (raw MSW) but constant discounting
# sccs.const
# sccs.ramsey

scc[5,'Mean'] =     mean(sccs.const) 
scc[5,'Median'] = median(sccs.const)  


# Scenario 6: Stochastic growth with Ramsey discounting
# incorporating the correlation
scc[6,'Mean'] =     mean(sccs.ramsey)  
scc[6,'Median'] = median(sccs.ramsey)  

# As a result, 95\% of all SCC draws fall between \$20 and \$220 per ton when using our 
# discounting rule, compared to the much wider range of \$15 to \$645 when using constant discounting.
# 95% confidence intervals for DICE
round(quantile(sccs.ramsey, probs=c(0.025, 0.05, 0.5, 0.95, 0.975)), 0)
round(quantile(sccs.const, probs=c(0.025, 0.05, 0.5, 0.95, 0.975)), 0)

scc

round(scc, 0)

scc.out = round(scc, 0)

rownames(scc.out) = c('1) IWG (DICE only, 2015 pulse, inflation adjusted)',
                      '2) 1 + Deterministic Emissions',
                      '3) 2 + Adjust Mean Population Growth to UW',
                      '4) 3 + Adjust Mean Economic Growth to MSW',
                      '5) 4 + Stochastic MSW Economic Growth',
                      paste0('6) 5 + Stochastic Discount Rate (rho=',round(rho*100,1),'%, eta=',round(eta,2),')'))

scc.out
library(scales)
usd <- dollar_format()
for (j in 1:ncol(scc.out)) scc.out[,j] = usd(as.numeric(scc.out[,j]))

library(stargazer)
stargazer(scc.out)

# ...all calibrated $\eta$ values for the 3\% rate shown in Table \ref{tab:results}.
# Incidentally, under all those $\eta$ values, the SCC estimates are all between \$60 and \$65
calc.scc(eta=output[['BR','3%']]$eta, rho=output[['BR','3%']]$rho, md=md.dice, g=md.g.dice)*mc.infl.adj
calc.scc(eta=output[['LVL','3%']]$eta, rho=output[['LVL','3%']]$rho, md=md.dice, g=md.g.dice)*mc.infl.adj
calc.scc(eta=output[['RW','3%']]$eta, rho=output[['RW','3%']]$rho, md=md.dice, g=md.g.dice)*mc.infl.adj
calc.scc(eta=output[['MR','3%']]$eta, rho=output[['MR','3%']]$rho, md=md.dice, g=md.g.dice)*mc.infl.adj
calc.scc(eta=output[['SS','3%']]$eta, rho=output[['SS','3%']]$rho, md=md.dice, g=md.g.dice)*mc.infl.adj
calc.scc(eta=output[['AADL','3%']]$eta, rho=output[['AADL','3%']]$rho, md=md.dice, g=md.g.dice)*mc.infl.adj


#### Figure 5. Plot 10k SCC values, under alternative discounting procedures (constant, Ramsey) -----
mean.g.dice = apply(md.g.dice, 2, weighted.mean, w=n.years.dice) # for each scenario, calculate LR average growth

dev.off()
# Plot 10k draws under DICE
my.cols = c('#ff6663','#50b161')
my.cols.alpha = scales::alpha(my.cols, 0.6)
my.cols[2] = colorspace::darken(my.cols[2], 0.5) # darken to differentiate in black and white
my.cols.alpha[2] = colorspace::darken(my.cols.alpha[2], 0.5) # darken to differentiate in black and white
sccs.const.plot = sccs.const
sccs.ramsey.plot = sccs.ramsey
library(scales)

pdf(figure.dir%&%'/SCC_MC_Draws_DICE_'%&%today()%&%'.pdf', family='serif', width=8.5, height=7)#, colormodel = 'gray')
par(mai=c(1.02, 0.82, 0.82, 1.5))
plot(log10(sccs.ramsey.plot) ~ I(100*mean.g.dice), type='n', 
     ylim=c(-2, 5), xlim=c(-2,5),
     cex.axis=1.2, cex.lab=1.5,
     main='', xlab='Time-averaged Growth Rate (%)', ylab='SCC ($/ton)', yaxt='n')
grid(nx=NULL, ny=NA)
ypts = 10^(-2:9)
axis(side=2, at=log10(ypts), labels=ypts, cex.axis=1.25)
abline(h=-2, v=0) # convention: horizontal axis at -2 in log10-space
abline(h=log10(ypts), col='lightgray', lty=3)
points(log10(sccs.const.plot) ~ I(100*mean.g.dice), col=my.cols.alpha[1], pch=1)
points(log10(sccs.ramsey.plot) ~ I(100*mean.g.dice), col=my.cols.alpha[2], pch=1)
abline(h=log10(mean(sccs.const.plot)), col=my.cols[1], lwd=2)
abline(h=log10(mean(sccs.ramsey.plot)), col=my.cols[2], lwd=2)
par(xpd=T)
text(x=5.3, y=log10(mean(sccs.ramsey.plot)), col=rgb(0.25,0.55,0.25), labels='E[SCC] = '%&%dollar(round(mean(sccs.ramsey.plot), 0)), pos=4, font=2)
text(x=5.3, y=log10(mean(sccs.const.plot)), col=my.cols[1], labels='E[SCC] = '%&%dollar(round(mean(sccs.const.plot), 0)), pos=4, font=2)
par(xpd=F)

legend(x=0, y=log10(10^-1), cex=1.2, pt.lwd=2, bty='n',
       legend=c('Constant 3% Discount Rate', 'Stochastic Discount Rate (With Correlation)'), 
       pch=1, col=my.cols)
rug(I(100*mean.g.dice), col=rgb(0,0,0,0.25), side=1)
rug(log10(sccs.const.plot), col=alpha(my.cols[1],0.25), side=2)
rug(log10(sccs.ramsey.plot), col=alpha(my.cols[2],0.25), side=2)
dev.off()





##### Run SCC calcs for a grid of eta values, varying rho to hold near-term r constant.
eta.grid = seq(0,2.1, by=0.025)
rho.grid = sapply(eta.grid, findrho, r.target=0.03) 

scc.eta.grid = mapply(FUN=function(e, r) calc.scc(eta=e, rho=r, md=md.dice, g=md.g.dice), eta.grid, rho.grid)*mc.infl.adj

scc.eta.grid[which(eta.grid==1.55)]
# for all $\eta$ values between 0.35 and 1.53, the discounting rule produces a mean SCC 
# that is between the deterministic SCC of \$56 (as in row 4) and the stochastic SCC of \$78 (in row 6).
# ...
# For example, the mean SCC exceeds \$100 for all $\eta$ values smaller than 0.12 or larger
# than 1.85 (again, varying $\rho$ to maintain a near-term 3\% discount rate).
tab = cbind(scc.eta.grid, eta.grid, rho.grid)

tab[between(tab[,1], scc[4,'Mean'], scc[6,'Mean']),] # eta values between $56 and $78

eta.rng = range(sapply(output[,'3%'], function(x) x$eta))
tab[between(tab[,2], eta.rng[1], eta.rng[2]),] # etas in the range of calibrated values for 3% rate (0.68 to 1.53)


mean(calc.scc(rho=findrho(1.53, r.target=0.03), eta=1.53, md=md.dice, g=md.g.dice, all=TRUE)*mc.infl.adj)
mean(calc.scc(rho=findrho(1.54, r.target=0.03), eta=1.54, md=md.dice, g=md.g.dice, all=TRUE)*mc.infl.adj)
tab[tab[,1]>100,]
mean(calc.scc(rho=findrho(0.13, r.target=0.03), eta=0.13, md=md.dice, g=md.g.dice, all=TRUE)*mc.infl.adj)
mean(calc.scc(rho=findrho(1.84, r.target=0.03), eta=1.84, md=md.dice, g=md.g.dice, all=TRUE)*mc.infl.adj)


calc.scc(eta=eta, rho=rho, md=md.dice, g=md.g.dice)*mc.infl.adj
calc.scc(eta=0, rho=0.03, md=md.dice, g=md.g.dice)*mc.infl.adj

eta.grid[which.min(scc.eta.grid)]

# Figure 4.
par(family='serif')
pdf(figure.dir%&%'/DICE_SCC_vs_eta_'%&%today()%&%'.pdf',
    family='serif', width=7, height=7)
plot(scc.eta.grid ~ eta.grid, type='l', lwd=2,  xlab=substitute(eta), ylab='SCC ($/ton)',
     ylim=c(0,200), main='', cex.main=1.5, cex.axis=1.5, cex.lab=1.5)
grid(col='gray70')
abline(h=scc[4,'Mean'], col=3, lwd=2, lty=3)
points(x=0, y=scc[4,1], pch=19, col=3, cex=1.5)
points(x=eta, y=scc[6,1], pch=19, col='black', cex=1.5)
legend('top', legend=c(as.expression(bquote('E[SCC], varying'~eta~'(and '*rho*')')),
                       as.expression(bquote('E[SCC] at '*eta==1.53*','~rho==0.8*'%')),
                       as.expression(bquote('Deterministic SCC,'~eta==0*','~rho==3*'%'))), lty=c(1,0,0), pch=c(NA,19,19), col=c('black','black',3),
       lwd=2, cex=1.3, bty='n')
dev.off()

